home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / b / b.lha / B / src / bint / b1nuQ.c < prev    next >
C/C++ Source or Header  |  1988-11-24  |  5KB  |  237 lines
  1. /* Copyright (c) Stichting Mathematisch Centrum, Amsterdam, 1985. */
  2.  
  3. /*
  4.   $Header: b1nuQ.c,v 1.4 85/08/22 16:51:40 timo Exp $
  5. */
  6.  
  7. #include "b.h"
  8. #include "b1obj.h"
  9. #include "b0con.h"
  10. #include "b1num.h"
  11.  
  12.  
  13. /* Product of integer and single "digit" */
  14.  
  15. Visible integer int1mul(v, n1) integer v; digit n1; {
  16.     integer a;
  17.     digit save, bigcarry, carry = 0;
  18.     double z, zz, n = n1;
  19.     register int i;
  20.     struct integer vv;
  21.  
  22.     FreezeSmallInt(v, vv);
  23.  
  24.     a = (integer) grab_num(Length(v)+2);
  25.  
  26.     for (i = 0; i < Length(v); ++i) {
  27.         z = Digit(v,i) * n;
  28.         bigcarry = zz = floor(z/BASE);
  29.         carry += z - zz*BASE;
  30.         Digit(a,i) = save = Modulo(carry, BASE);
  31.         carry = (carry-save)/BASE + bigcarry;
  32.     }
  33.  
  34.     Digit(a,i) = save = Modulo(carry, BASE);
  35.     Digit(a,i+1) = (carry-save)/BASE;
  36.  
  37.     return int_canon(a);
  38. }
  39.  
  40.  
  41. /* Quotient of positive integer and single "digit" > 0 */
  42.  
  43. Hidden integer int1div(v, n1, prem) integer v; digit n1, *prem; {
  44.     integer q;
  45.     double r_over_n, r = 0, n = n1;
  46.     register int i;
  47.     struct integer vv;
  48.  
  49.     FreezeSmallInt(v, vv);
  50.  
  51.     q = (integer) grab_num(Length(v));
  52.     for (i = Length(v)-1; i >= 0; --i) {
  53.         r = r*BASE + Digit(v,i);
  54.         Digit(q,i) = r_over_n = floor(r/n);
  55.         r -= r_over_n * n;
  56.     }
  57.     if (prem)
  58.         *prem = r;
  59.     return int_canon(q);
  60. }
  61.  
  62.  
  63. /* Long division routine, gives access to division algorithm. */
  64.  
  65. Visible digit int_ldiv(v1, w1, pquot, prem) integer v1, w1, *pquot, *prem; {
  66.     integer a;
  67.     int sign = 1, rel_v = 0, rel_w = 0;
  68.     digit div, rem;
  69.     struct integer vv1, ww1;
  70.  
  71.     if (w1 == int_0) syserr(MESS(1100, "zero division (int_ldiv)"));
  72.  
  73.     /* Make v, w positive */
  74.     if (Msd(v1) < 0) {
  75.         sign = -1;
  76.         ++rel_v;
  77.         v1 = int_neg(v1);
  78.     }
  79.  
  80.     if (Msd(w1) < 0) {
  81.         sign *= -1;
  82.         ++rel_w;
  83.         w1 = int_neg(w1);
  84.     }
  85.     
  86.     FreezeSmallInt(v1, vv1);
  87.     FreezeSmallInt(w1, ww1);
  88.  
  89.     div = sign;
  90.  
  91.     /* Check v << w or single-digit w */
  92.     if (Length(v1) < Length(w1)
  93.         || Length(v1) == Length(w1)
  94.             && Digit(v1, Length(v1)-1) < Digit(w1, Length(w1)-1)) {
  95.         a = int_0;
  96.         if (prem) {
  97.             if (v1 == &vv1) *prem= (integer) MkSmallInt(Digit(v1,0));
  98.             else *prem = (integer) Copy(v1);
  99.         }
  100.     }
  101.     else if (Length(w1) == 1) {
  102.         /* Single-precision division */
  103.         a = int1div(v1, Digit(w1,0), &rem);
  104.         if (prem) *prem = mk_int((double)rem);
  105.     }
  106.     else {
  107.         /* Multi-precision division */
  108.         /* Cf. Knuth II Sec. 4.3.1. Algorithm D */
  109.         /* Note that we count in the reverse direction (not easier!) */
  110.  
  111.         double z, zz;
  112.         digit carry, save, bigcarry;
  113.         double q, d = BASE/(Digit(w1, Length(w1)-1)+1);
  114.         register int i, j, k;
  115.         integer v, w;
  116.         digit vj;
  117.  
  118.         /* Normalize: make Msd(w) >= BASE/2 by multiplying
  119.            both v and w by d */
  120.  
  121.         v = int1mul(v1, (digit)d);
  122.             /* v is used as accumulator, must make a copy */
  123.             /* v cannot be int_1 */
  124.             /* (then it would be one of the cases above) */
  125.  
  126.         if (d == 1) w = (integer) Copy(w1);
  127.         else w = int1mul(w1, (digit)d);
  128.  
  129.         a = (integer) grab_num(Length(v1)-Length(w)+1);
  130.  
  131.         /* Division loop */
  132.  
  133.         for (j = Length(v1), k = Length(a)-1; k >= 0; --j, --k) {
  134.             vj = j >= Length(v) ? 0 : Digit(v,j);
  135.  
  136.             /* Find trial digit */
  137.  
  138.             if (vj == Digit(w, Length(w)-1)) q = BASE-1;
  139.             else q = floor( ((double)vj*BASE
  140.                     + Digit(v,j-1)) / Digit(w, Length(w)-1) );
  141.  
  142.             /* Correct trial digit */
  143.  
  144.             while (Digit(w,Length(w)-2) * q >
  145.                 ((double)vj*BASE + Digit(v,j-1)
  146.                     - q*Digit(w, Length(w)-1)) *BASE + Digit(v,j-2))
  147.                 --q;
  148.  
  149.             /* Subtract q*w from v */
  150.  
  151.             carry = 0;
  152.             for (i = 0; i < Length(w) && i+k < Length(v); ++i) {
  153.                 z = Digit(w,i) * q;
  154.                 bigcarry = zz = floor(z/BASE);
  155.                 carry += Digit(v,i+k) - z + zz*BASE;
  156.                 Digit(v,i+k) =
  157.                     save = Modulo(carry, BASE);
  158.                 carry = (carry-save)/BASE - bigcarry;
  159.             }
  160.  
  161.             if (i+k < Length(v))
  162.                 carry += Digit(v, i+k), Digit(v, i+k) = 0;
  163.  
  164.             /* Add back necessary? */
  165.  
  166.                 /* It is very difficult to find test cases
  167.                    where add back is necessary if BASE is large.
  168.                    Thanks to Arjen Lenstra, we have v=n*n-1, w=n,
  169.                    where n = 8109636009903000000 (the last six
  170.                    digits are not important). */
  171.  
  172.             if (carry == 0)        /* No */
  173.                 Digit(a,k) = q;
  174.             else {        /* Yes, add back */
  175.                 if (carry != -1) syserr(MESS(1101, "int_ldiv internal failure"));
  176.                 Digit(a,k) = q-1;
  177.                 carry = 0;
  178.                 for (i = 0; i < Length(w) && i+k < Length(v); ++i) {
  179.                     carry += Digit(v, i+k) + Digit(w,i);
  180.                     Digit(v,i+k) =
  181.                         save = Modulo(carry, BASE);
  182.                     carry = (carry-save)/BASE;
  183.                 }
  184.             }
  185.         }    /* End for(j) */
  186.  
  187.         if (prem) *prem = int_canon(v);    /* Store remainder */
  188.         else release((value) v);
  189.         div = sign*d;    /* Store normalization factor */
  190.         release((value) w);
  191.         a = int_canon(a);
  192.     }
  193.  
  194.     if (rel_v) release((value) v1);
  195.     if (rel_w) release((value) w1);
  196.  
  197.     if (sign < 0) {
  198.         integer temp = a;
  199.         a = int_neg(a);
  200.         release((value) temp);
  201.     }
  202.  
  203.     if (pquot) *pquot = a;
  204.     else release((value) a);
  205.     return div;
  206. }
  207.  
  208.  
  209. Visible integer int_quot(v, w) integer v, w; {
  210.     integer quo;
  211.     VOID int_ldiv(v, w, &quo, (integer*)0);
  212.     return quo;
  213. }
  214.  
  215. Visible integer int_mod(v, w) integer v, w; {
  216.     integer rem;
  217.     digit div;
  218.     bool flag;
  219.     div = int_ldiv(v, w, (integer*)0, &rem); /* Rem. is always positive */
  220.     if (rem == int_0)
  221.         return rem; /* v mod w = 0 */
  222.     flag = (div < 0);
  223.     if (flag || Msd(w) < 0) div = -div;
  224.     if (div != 1) {    /* Divide by div to get proper remainder back */
  225.         v = int1div(rem, div, (digit*)0);
  226.         release((value) rem);
  227.         rem = v;
  228.     }
  229.     if (flag) { /* Make same sign as w */
  230.         if (Msd(w) < 0) v = int_sum(w, rem);
  231.         else v = int_diff(w, rem);
  232.         release((value) rem);
  233.         rem = v;
  234.     }
  235.     return rem;
  236. }
  237.